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Abstract 

We study effects of nonlocality of the cubic self-focusing nonlinearity on the stability and 
symmetry-breaking bifurcation (SBB) of solitons in the model of a planar dual-core optical waveg- 
uide with nonlocal (thermal) nonlinearity. In comparison with the well-known coupled systems 
with the local nonlinearity, the present setting is affected by the competition of different spatial 
scales, viz., the coupling length and correlation radius of the nonlocality, \fd. By means of numer- 
ical methods and variational approximation (VA, which is relevant for small d), we find that, with 
the increase of the correlation radius, the SBB changes from subcritical into supercritical, which 
makes all the asymmetric solitons stable. On the other hand, the nonlocality has little influence on 
the stability of antisymmetric solitons. Analytical results for the SBB are also obtained (actually, 
for antisymmetric "accessible solitons") in the opposite limit of the ultra- nonlocal nonlinearity, 
using a coupler based on the Snyder-Mitchell model. The results help to grasp the general picture 
of the symmetry breaking in nonlocal couplers. 

PACS numbers: 42.65.Tg, 42.65.Jx, 42.65.Wi, 03.75.Lm 
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I. INTRODUCTION 



Dual-core systems, featuring intrinsic nonlinearity in parallel cores coupled by linear tun- 
neling of wave fields, find their realizations in various physical settings. Well-known systems 
of this type in optics are twin-core fibers (see also an early review |7( ) and Bragg grat- 

ings [8], as well as double planar waveguides with the second-harmonic-generating intrinsic 
nonlinearity 
condensates 



9|. Similar settings for matter waves are represented by two-layer Bose-Einstein 



10|- 42|. A fundamental physical effect in nonlinear symmetric dual-core sys- 



tems is the symmetry-breaking bifurcation (SBB), alias the phase transition, which desta- 
bilizes symmetric modes and gives rise to asymmetric ones. In nonlinear optics, the SBB 
was studied in detail for continuous- wave (spatially uniform) states [6| and solitons [8j, [131 in 

nn n □ n 

twin-core fibers |4|, |5], [12]-[17[ and Bragg gratings [8] with the Kerr (cubic) nonlinearity, as 



well as for solitons in double-core waveguides with the quadratic [9| and cubic-quintic [18 1 



non 
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inearity. The SBB was studied too for matter-wave solitons in two-layer BEC settings 



111- 



The self-focusing cubic nonlinearity gives rise to the SBB of the subcritical type (alias 
the phase transition of the first kind) for solitons in the symmetric dual-core system. The 
bifurcation of this type is characterized by originally unstable branches of emerging asym- 
metric modes, which at first extend backward (in the direction of weaker nonlinearity), and 



then turn forward, retrieving the stability at the turning points [19| . In this case, the system 
demonstrates a bistability and hysteresis in a limited interval, characteristic to phase transi- 
tions of the first kind. If the dual-core system is equipped with a periodic potential (lattice) 
acting in the direction transversal to the propagation coordinate, the character of the SBB 



changes to supercritical above a certain threshold value of the lattice's strength 



ll|. The 



supercritical bifurcation (alias the phase transition of the second kind) gives rise to stable 
branches of asymmetric modes going in the forward direction 19fl. The SBBs belong to this 

Tin 

type too in the twin-core Bragg grating, and in quadratically nonlinear waveguides [8|, |9J . 
In addition to numerical analysis of symmetric, antisymmetric and asymmetric soliton 
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modes in dual-core system with the intrinsic cubic nonlinearity 13|, |l5), the bifurcation point 



was found in an exact analytical form 



4J , and the emerging asymmetric solitons were studied 



12 



14, hj. The latter method 



in detail by means of the variational approximation (VA) 
is relevant for studies of solitons in many models originating in nonlinear optics and related 
fields 20], while the possibility to find the exact bifurcation point is a feature specific to 



particular systems. 

The nonlinear response in optical media may feature spatial nonlocality, which means 



that the local change of the refractive index induced by the 



distribution of the light intensity in a vicinity of a given point 



21 



ight beam depends on the 



221 ] . The nonlocality arises 



when the nonlinear optical response involves mechanisms such as heat diffusion, as analyzed 



theoretically 
crystals |26|, 



23J and demonstrated experimentally 24j,|25|], molecular reorientation in liquid 



271 ]. atomic diffusion 28M30| . etc. The fields of nanophotonics and plasmonics 



also give rise to effective nonlocalities, due to light-matter interactions occurring in these 
media on deeply subwavelength scales [3l|, [32]. 

Nonlocal nonlinearities are known in other physical media, including plasmas 33J and 
self-gravitating photonic beams 34j. Long-range interactions play an important role in 
dipolar Bose-Einstein condensates (BECs) too 35], and nonlocal gravity-like interactions 
can be induced in BEC by means of laser illumination 36 1. 

The nonlocality, which introduces a new spatial scale, namely, the correlation radius 
(denoted below as \/d), may drastically alter nonlinear excitations in optical systems, due 
to the interplay of y/dwiih other natural scales. In particular, the nonlocality changes 
the character of interactions between solitons 37], and it suppresses the beam's collapse and 



transverse instabi 
of soliton modes 



ities 



38 



41 



391 ] . The nonlocality also accounts for the formation of new types 



However, to the best of knowledge, the influence of nonlocality 
on the performance of optical couplers has not been reported yet. In particular, new effects 
may be expected due to the competition of \fd with the coupling length, i.e., the interplay 
of nonlocal and local interactions. This is the objective of the present work. 

We consider the formation of solitons in a planar dual-core waveguide, in which the 
nonlocal nonlinearity of the thermal type acts in both cores, while the coupling between 
them remains linear and local, as the heat diffusion does not transfer energy across the 
gap separating the waveguides. Similar to couplers with the local nonlinearity, the nonlocal 
model gives rise to three types of solitons, viz., symmetric, antisymmetric and asymmetric 
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ones. However, the nonlocality significantly affects the symmetry-breaking phase transition 
(SBB) for solitons, as well as stability of the emerging asymmetric solitons, which are basic 
properties of nonlinear couplers: at a critical value of the \fd, the SBB changes its character 
from sub- to supercritical. Taking into regard the potential that nonlinear couplers have for 
various application to photonics, such as all-optical switching {2, 2|, the use of the nonlocality 
for the control of the soliton dynamics in these systems may help to expand the range of 
the applications. While our analysis is performed in terms of the thermal nonlinearity in 
optical waveguides, the results may plausibly apply to other dual-core physical systems 
which feature the nonlocal nonlinearity. 

The paper is organized as follows. The model is formulated in Section II, and analytical 
results are reported in Section III. These results are obtained by means of the VA for 
solitons in the case of weak nonlocality (small y/d), and, on the other hand, the SBB is also 
investigated (in fact, for antisymmetric solitons) in the opposite limit of the ultra-nonlocal 
nonlinearity, in terms of a coupled system for "accessible solitons" [the Snyder-Mitchell (SM) 



model 2l|]]. In particular, the exact bifurcation point is found for the SM system. The results 
for the small correlation radius explicitly demonstrate the shift of the SBB point to larger 
values of the soliton's power, and the trend to the transition of the subcritical bifurcation 
into the supercritical one, while the findings reported for the ultra-nonlocal system help 
to apprehend the general situation. Numerical results, which provide the full description 
of solitons in the nonlocal dual-core system for moderate values of the correlation radius, 
are presented in Section IV. In the case of the weak nonlocality, these results verify the 
analytical results produced by the VA. The paper is concluded by Section V. 

II. THE MODEL 

The propagation of optical beams along axis z in the planar dual-core waveguide with 



the intrinsic self-focusing nonlinearity of the thermal type 



22 



23| is described by the system 



of linearly coupled nonlinear Schrodinger (NLS) equations for complex field amplitudes u, v 
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in the two cores, and respective local perturbations m, n of the refractive index: 

iu z + ^u xx + mu + v — 0, (la) 

iv z + ^v xx + nv + u = 0, (lb) 

m — dm xx = \u\ 2 , (lc) 
n - dn xx = \v\ 2 , (Id) 

where x is the transverse coordinate, the coupling constant [the coefficient in front of terms 
v and u in Eqs. (pji) and ([lb), respectively] is scaled to be 1 (accordingly, the coupling 
length is also ~ 1), and d is the squared correlation radius of the nonlocality. In fact, d 
controls the competition between the length scales determined by the nonlocal and local 
interactions in the system. 

Stationary solutions to Eqs. (0Q) with propagation constant b are looked for as 

u (z, x) = e ibz U(x),v (z, x) = e ibz V(x), (2a) 
m = m(x),n = n(x), (2b) 

with real functions U(x) and V(x) obeying the following equations: 

-bU + -U" + mU + V = 0, (3a) 

-bV + -V" + nV + U = 0, (3b) 

m - dm" = U 2 , (3c) 

n - dn" = V 2 . (3d) 

Equations (IT]) conserve the total power, 

/oo roc 
\u\ 2 dx+ / \v\ 2 dx. (4) 
-oo J — oo 

Obviously, symmetric [U(x) = V(x)] and antisymmetric [U(x) = —V(x)] modes have Pi = 
P2, while asymmetric ones can be characterized by parameter 

which takes values — 1 < 6 < +1. 
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Parallel to Eqs. (JTJ) , it is relevant to consider the ultra-nonlocal model, taken in the form 



of two linearly coupled SM equations 



2l|, 



1 1 2 

i u z + -u xx --P u xu + v = 0, 



iv z + -=v x 



1 



-P,,x v + u = 0, 



(6a) 
(6b) 



where P\p are the powers defined as per Eq. (@|. Actually, Eqs. (J6j) correspond to the 
version of Eqs. ([lb) and ((Tbl) with spatially averaged right-hand sides. To the best of 
our knowledge, the SM coupler was not considered before, while the extreme nonlocality 
postulated in the SM model per se finds realizations and applications in diverse optical [43j 
and optomechanical 44] settings. 



III. ANALYTICAL RESULTS 



A. The variational approximation for the weakly nonlocal system 



To apply the VA to the present system, we note that, in the case of weak nonlocality 
(d <C 1), Eqs. ([3b) and (EU) yield, in the first approximation, m = U 2 + d(U 2 )", and 
n = V 2 + d (V 2 )" [45]. The substitution of this approximation into Eqs. and ([3b) leads 
to a system of two coupled equations with nonlinear-diffraction terms: 



-bU + hj" + U 3 + dU (U 2 )" + V = 0, 
-bV + \v" + V 3 + dV (V 2 )" + U = 0, 



which may be derived from the Lagrangian with density 
C 



1 

4 
+d 



(Uf + [V'f] + \ (U 2 + V 2 ) -\{u* + V) 



U 2 {U'f + V 2 {V'f 



uv. 



The ansatz for soliton solutions may be naturally chosen as 



{U(x), V{x)} = {A, B} sech (x/W) 



(7a) 
(7b) 



(8) 



(9) 



where A and B are amplitudes of the two components, and W is their common width. 
The substitution of the ansatz into density (JHJ) and evaluation of the integrals yields the 
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corresponding Lagrangian, 

/+oo A2 , r>2 

^ Cdx = ^^ + b(A 2 + B 2 )W 

-I (A 4 + P 4 ) W + Ad{ ^lt B4) - 2ABW. (10) 

This Lagrangian can be more conveniently rewritten in terms of the total power P, see Eq. 
flU), and power imbalance Q = Pi — P 2 , 

2 (A 2 + B 2 )W = P, 2 (A 2 -B 2 )W = Q, (11) 

as follows: 

p P 2 + n 2 

4f^f-o^?, (13) 

where a = 1 for symmetric solitons and asymmetric ones generated from them by the SBB, 
and a = — 1 for antisymmetric solitons. The corresponding Euler-Lagrange equations are 
OL/dW = dL/dQ = OL/dP = 0, i.e., 

P P 2 + Q 2 3dP 2 + Q 2 



W 4 5 W 2 



0, (14) 



Q {-^ + ^ + ^w^wj=^ (15) 

_P 1 2dP aP 

~QW~ QW 2 ~ 15W 3 + ^/P2 _ Q2 ( > 

Equation ffl6|) . which determines the propagation constant, b, is detached from Eqs. 
( fli"]) and (fT5l) . Equation ( Tl~5l) yields either Q = 0, which corresponds to symmetric and 
antisymmetric solitons, or 

1 2d 1 



15W 7 " 3 ^P 2 -Q 

for asymmetric ones. Further, the expansion of Eqs. ( JT4l and ( fT6l) for small <i, i.e., the weak 
nonlocality, yields 

W ^ — + —P, b « — P 2 + a - — P 4 , (18) 
P 5 ' 32 192 ' K J 

which predicts that, naturally, the nonlocality makes the soliton wider, for given total power 
P. This is confirmed by the numerical solutions, as shown below. 



The most essential point is to find the critical power Pbif, at which the asymmetric solitons 
bifurcate from the symmetric ones. This value is determined by a system of equations ( I14p 
and (1171) . in which one should set Q — 0. Further, using the assumption of the weak 
nonlocality, i.e., small d, the ensuing solution for P^f can be expanded up to order d, which 
yields 

P bif = 2VQ+ (24V6/5) d. (19) 



Note that, at d = 0, Eq. (j!9p gives P m (d = 0) = 2^6 [12], which may be compared to the 
known exact result [4j, (Pbif) exact = 8/v3, the relative error being 0.057. 

The VA predicts, as per Eq. ( Fl9l) . the increase of the soliton's power at the bifurcation 
point due to the weak nonlocality. To compare the prediction with the numerical findings, 
we take the slope of the Pbif(^) dependence at d — 0, for which Eq. ffT9]) yields 

~d(P U{ ) 



\d=0 



= 24^6/5 « 11.758. (20) 

variational 



d(d) 

On the other hand, the same slope obtained from the numerical solution (see the next 
section) is 

^§r ld=o l ^ 10 - 867 ' (21) 

\ J J numerical 

the relative error of the VA prediction being 0.075 (see Table 1). 

It is also possible to find another critical power, P t h, which corresponds to the turn- 
ing point (i.e., the stabilization threshold for asymmetric solitons) on the dependence of 
the asymmetry parameter, O = Q/P [see Eq. (j5J)], on total power P. To this end, one 
should obtain a dependence between and P, eliminating W from Eqs. f TH|) and ffl7]) . and 
identifying P t h from condition 



In the limit of d = 0, the result produced by the VA is known [171 ]: 



(Pth), =0 = 3 ■ 6 1/4 ^ 4.695, (23) 

the corresponding value of the asymmetry at the critical point being Oth = 1/ v3- On the 
other hand, the numerically found threshold power at d = is 

[(^rXum * 4.548, (24) 

hence the relative error produced by the comparison of Eqs. (|23|) and (I2"4"j) is 0.031 (see 
Table 1). 
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TABLE I The comparison between the VA-predicted characteristics of the symmetry-breaking 
bifurcation, in the local and weakly nonlocal systems, and their numerically found counterparts. 



Parameter 


VA 


Numeric 


VA-Numer /o/ \ 
VA ( /o ) 


-Pth|d=0 


4.6953 


4.5484 


3.13 


Pbif\d=0 


4.8989 


4.6188 


5.72 


d{P th ) | 
d(d) l d =o 


12.2681 


13.8270 


-12.71 


d(P bii ) | 
d(d) \d=0 


11.7575 


10.8666 


7.58 



Further, the expansion of Eqs. (1141) . ( fl~Tj) and ( 1221) for small d yields the following pre- 
diction for the slope of curve Pth(d) at d = 0: 



d(Pt 



th 



d(d) 



\d=0 



16 



6 3/4 « 12.268, 



(25) 



variational 



while the numerically found counterpart of this value is 

~d(P th ) . 



_ d(d) 



\d=0 



13.827, 



(26) 



hence the respective relative error is 0.127 (see Table 1). 
Finally, we note that the relation 

^r |d=o> ^r |d=o ' (27) 

see Eqs. fl25|) and fl20|) . suggests that P t h and Pbif will eventually merge into a single 
critical/threshold value, which implies the transition from the subcritical bifurcation to the 
supercritical one, as confirmed by numerical results displayed below. 



B. The coupler for "accessible solitons" (the Snyder-Mitchell model) 



In the opposite case of the ultra-nonlocal nonlinearity, substitution (T5k) transforms cou- 
pled SM equations 21] and Eq. (BJ into their stationary versions: 



-bU + \u" - \p u x 2 \J + V 



0. 



-bV + \v" - \p v x 2 V + U = 0, 



(28a) 
(28b) 
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/OO POO 
U 2 (x)dx, P v = V 2 (x)dx. (29) 
-oo J — OO 

In spite of the apparently simple form of Eqs. ( 1281) and ([29]) . it is not possible to find 
exact solutions for asymmetric solitons. A solution can be obtained, by means of the WKB 
approximation, in the limit case of the strong asymmetry, P v <C P u . In this case, the 
U component is tantamount to the ground-state wave function of the harmonic oscillator 
(HO), with the corresponding HO length L u = Pu 1 ^, eigenvalue of the propagation constant 
b = —sj P u /2, and amplitude 

U(x = 0) = tt- 1 /^ 8 , (30) 

while the weak V component develops a broad shape, with a small amplitude, V(x = 
0) ~ —a/2/ ttPv P u 1 ^ 8 , and large width, L v 2a/ y/P^/P v . The wave function of the V- 
component can be written in a relatively simple explicit WKB form in the "resonant" case, 



P V = P U /(2(2N+1)Y, (31) 

with large integer N, when the (2iV)-th energy eigenvalue in the shallow HO potential 
(assuming that N = corresponds to the ground state) in the ^-component is matched to 
the ground-state eigenvalue of the HO in the £7-component: 

1/4 



V{x) 



x cos < 



fp, 



■ arcsm 




(32) 



at x 2 < yfP^/P V) and V(x) = at x 2 > y/P^/P v [if resonance condition ( 13T]) does not hold, 
the WKB expression fl32|) needs a correction around the edge points, x 2 = y/P^/P v ]. 

It follows from Eq. (J28b) taken at the inflexion point (V" = 0) closest to x = that 
the strongly asymmetric mode has opposite signs of U(x = 0) and V(x = 0) (as written in 
the above formulas), i.e., this asymmetric state develops from the antisymmetric one. The 
respective point of the antisymmetry-breaking bifurcation can be found in an exact form. 
To this end, a solution to Eqs. ( 128]) near the bifurcation point is looked for as 



{U(x),V(x)} = ±U exp[ -\jjx 2 ) +5U(x) 



(33) 
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where the propagation constant and amplitude of the lowest 
mode, with total power P (in both components), are 

6=-l-(l/2)/P/2, 



unperturbed 



antisymmetric 
(34) 



U = vr" 1 / 4 (P/2) 5/8 , (35) 

cf. Eq. (!30|) . and an infinitesimal antisymmetry-breaking perturbation, 5U(x), obeys the 
equation following from the substitution of expression fl33l into Eqs. ([28]) and ff29|) and 
subsequent linearization: 

(l-b)6U + ^6U" -^Px 2 5U = U (5P)x 2 exp [~\\[^ x ^\ > ( 36a ) 
5P = U I exp I --\l-x 2 \ 5U{x)dx. (36b) 




2 V 2 

A relevant solution to inhomogeneous equation (136k) can be found as 

5U = (d + U 2 x 2 ^j exp [~\^ 2 \ , (37a) 

U 6P 4U 5P 

On = ; , 02 = ; • (37b) 

3yfP/2-4 3yfp/2-4 V ; 

Finally, substituting expressions ( 157|) into Eq. (I3"6"b) and canceling o"P as a common factor, 
the self-consistency condition yields a simple exact result for the total power at which the 
increase of the spontaneous breaking of the antisymmetry occurs: p c ( r antls y mm ) = g # 

IV. NUMERICAL RESULTS 

Numerical solution of Eqs. ()3]) was performed by means of the standard relaxation 
method. As predicted by the VA, three soliton families, symmetric, asymmetric, and an- 
tisymmetric ones, persist in the nonlocal system. The numerically found relation between 
the total power, P, and propagation constant b for symmetric and antisymmetric solutions 
is shown in Figs. [TJ It is seen that b monotonically grows with P at a fixed value of the 
nonlocality range, \/o![which implies that the solitons may be stable in terms of the Vakhitov- 



Kolokolov (VK) criterion 



46)], and b decreases with d at fixed P. Both these properties 



are correctly predicted by the VA, see Eq. ffT8]) . The fact that all the curves originate, 
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at P = 0, from the same point, is obvious, as it immediately follows from Eqs. ([3]) that 
limp^o b(P) = a = sgn (UV). 




FIG. 1. (Color online) Total power P versus the soliton's propagation constant, b, at different 
fixed values of the squared nonlocality correlation radius, d, for symmetric solitons in the model 
based on Eqs. ([I]). The inset shows a typical soliton profile. For the antisymmetric solitons, b is 
shifted by Ab = — 2. All quantities are plotted in arbitrary dimensionless units. 



Proceeding to numerically found asymmetric solitons, in Fig. [21(a) we plot the respective 
P(b) curves for for different fixed values of d. As in the local system, asymmetric modes 
appear through the SBB when the total power exceeds the threshold value, P t h- Note that 
the threshold, as well as the value of the total power at the bifurcation point, P = P^, 
significantly grow with d [see Fig. 121(b)], in accordance with the prediction of the VA given 
by Eqs. ( |23i) and ( 120]) . Further, the P(b) curves change their shape with the growth of 
the nonlocality radius: At small d, the slope, dP/db, is initially negative (which definitely 
implies the instability, according to the VK criterion 46[), going over to dP/db > with the 
further increase of b. With the increase of d, the segment with the negative slope shrinks, 
and disappears at d > 0.05. 

The change in the shape of the P(b) characteristics is directly related to the switch of 
the SBB from sub- to the supercritical type (in other words, the switch from the symmetry- 
breaking phase transition from the first to second kind) [191 ]. as shown in Fig. 12(c), where 
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P = P t h determines the turning points of the 0(P) curves, and their unstable portions with 
dQ/dP < precisely correspond to the segments with dP/db < in Fig. 2(b), both being 
confined to P t h < P < Pbif- Accordingly, the type of the SBB is subcritical, with P th < P bif 
at d < 0.05, and supercritical, with P t h = P^f, at d > 0.05. The merger of P t h and P^f 
into the single value at d > 0.05 is clearly observed in Fig. [2]^b). Recall that, as mentioned 
above, the trend to the merger of the two critical powers was predicted by the VA, see Eq. 

(E7D- 

It is relevant to compare this result with the transition from the subcritical SBB for 
solitons into the supercritical bifurcation under the action of the periodic potential 11]. 
Although the models are very different (the one considered in Ref. ll| is local), a common 
feature is the introduction of a specific spatial scale — the nonlocality range in the present 
model, \/d, or the lattice period in the local model — which is a factor accounting for the 
change of the character of the SBB. 

The stability of the solitons was tested by means of systematic simulations of Eqs. (CD), 
starting with perturbed initial conditions, u(x, z = 0) = U(x)(l + p(x)),v(x, z = 0) = 
V(x)(l + p(x)), where U(x),V(x) is a stationary solution, and p(x) is a small-amplitude 
random function. As expected, it has been found that the solid portions of the curves in 
Figs. EJ^a) and[2^c), with dP/db > and dQ/dP > 0, carry stable solitons, while the dashed 
segments, with dP/db < and dQ/dP < 0, represent unstable solutions. Thus, the increase 
of the nonlocality radius, \/d, gradually eliminates the instability region for the asymmetric 
solitons, making them completely stable in the case when the SBB is supercritical, i.e., at 
d > 0.05. 

It is relevant to explore the evolution of the two species of unstable solitons in the dual- 
core system, viz., asymmetric ones belonging to the segments of the 0(P) curves with the 
negative slope [i.e., < 0(P t h), that exist at d < 0.05], which are represented, for example, 
by point B in Fig. |2]^c), and symmetric solitons with P > Pbif, sampled by point D in Fig. 
EJ^c). Figure G^a) displays the result for the unstable asymmetric soliton, which demonstrates 
long-lived oscillations, initiated by the instability, and eventual relaxation into a stable 
soliton with almost the same power but higher asymmetry, > 0(P t h), which belongs 
to the stable branch of asymmetric modes in Fig. EJ^c). Further, Fig. [3]^b) demonstrates 
that the instability of the symmetric soliton leads to its spontaneous rearrangement into an 
asymmetric one, with nearly the same total power. 
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FIG. 2. (Color online) (a) Total power P versus propagation constant b for asymmetric solitons 
at different values of the squared nonlocality radius, d. The inset shows a typical soliton profile. 

(b) The dependence on nonlocality d of the total power, Pbif ; at which the symmetry-breaking 
bifurcation gives rise to asymmetric solitons, and of the threshold power, P t h, at which the pair 
of stable and unstable asymmetric solitons emerge subcritically. (c) The bifurcation diagram 
accounting for the creation of the asymmetric solitons from the symmetric ones. In panels (a) and 

(c) , dashed curves depict unstable portions of the asymmetric-soliton families [the border between 
stable and unstable (dashed) parts of the symmetric-soliton family in Fig. 2(c) corresponds to 
d = 0.01]. All quantities are plotted in arbitrary dimensionless units. 

We have also studied the stability and evolution of antisymmetric solitons for different 
strengths of the nonlocality in the model based on Eqs. ([1]) (the stability of the antisymmetric 
solitons in the model of the coupler with the local nonlinearity was studied, in a numerical 
form, in Ref. [jjj]). In contrast to the asymmetric solitons, where the nonlocality leads to the 
transition from the subcritical SBB to the supercritical bifurcation, and thus enhances the 



14 



Oh 



< 

X 

03 



1.8 
1.6 H 
1.4 
1.2 



i — 1 — i — 1 — i — <—//- 
50 100150 



1.1 
0.9 
0.7 

0.6 
0.4 
0.2 
0.0 




i 1 1 1 1 1 r 

50 100150 



-//- 




1 — i — 1 — i — 1 — i — '—//- 
50 100150 



< 

x 

03 



1.2 -f 

0.8 

0.4 



-//- 







50 



Oh 
B 
< 

X 
03 



2.1- 
1.6- 
1.1 




-//- 







50 



1.0 
© 0.5 
0.0 




-7A 



50 



(a) 




1000 

z 



1000 

z 



1000 

z 



(b) 



1000 

z 



1000 

z 



1000 

z 



Max Amp(u) 



2000 3000 
Max Amp(v) 



2000 



3000 

-0 



2000 



3000 



Max Amp(u) 



2000 3000 
— Max Amp(v) 



2000 



3000 




2000 



3000 



FIG. 3. (Color online) (a,b): The evolution of perturbed unstable solitons corresponding, re- 
spectively, to points B and D marked in Fig. 2(c) (weakly asymmetric and symmetric solitons) is 
shown in terms of amplitudes of both components, and asymmetry measure ([5]). Both examples 
pertain to d = 0.01. All quantities are plotted in arbitrary dimensionless units. 
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stability of the asymmetric solitons, it has been found that the stability of the antisymmetric 
ones is weakly affected by the nonlocality: the stability region slightly expands under the 
action of the nonlocality, without dramatic changes. 

V. CONCLUSION 

We have introduced the nonlocal generalizations of the standard model of the nonlinear 
directional coupler. The system can be built, in particular, as a dual-core optical waveguide 
made of a material with thermal nonlinearity. By means of the VA (variational approxima- 
tion) and systematic numerical analysis, we have found that the relatively weak nonlocality 
shifts the SBB (symmetry-breaking bifurcation) of solitons to larger values of the total 
power, and eventually changes the character of the SBB from subcritical to the supercritical 
(i.e., the corresponding phase transition of the first kind goes over into the transition of 
the second kind). Thus, the nonlocality of the cubic nonlinearity enhances the stability for 
the asymmetric solitons, and eventually leads to their stabilization in the whole existence 
domain, while only slightly affecting the stability of antisymmetric solitons. For the con- 
sideration of the opposite case of the ultra-nonlocal nonlinearity, the coupler based on the 
SM (Snyder-Mitchell) model was introduced. In that case, the phase transition leads to the 
spontaneous breaking of the antisymmetry of the corresponding two-component "accessible 
solitons". The exact transition point was found, and the strongly asymmetric states were 
found by means of the WKB approximation. 

The analysis reported in this paper can be extended in other directions. In particular, as 
concerns nonlocal dual-core systems in other physical contexts, it may be quite interesting 
to study the SBB and asymmetric solitons in the case when the nonlocal interactions act 
between the cores, an important example being a two-layer dipolar BEC The symmetry- 
breaking point can be easily found for the respectively modified SM coupler model. A 
challenging extension is to construct two-dimensional solitons in dual-core systems, where 
they may be stabilized against the collapse by the nonlocality of the nonlinearity. 
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